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ABSTRACT 

Recent numerical simulations have shown long-lived axisymmetric sub- and super-Keplerian flows 
in protoplanetary disks. These zonal flows are found in local as well as global simulations of disks 
unstable to the magnetorotational instability. This paper covers our study of the strength and lifetime 
of zonal flows and the resulting long-lived gas over- and underdensities as functions of the azimuthal 
and radial size of the local shearing box. We further investigate dust particle concentrations without 
feedback on the gas and without self-gravity. Strength and lifetime of zonal flows increases with 
the radial extent of the simulation box, but decreases with the azimuthal box size. Our simulations 
support earlier results that zonal flows have a natural radial length scale of 5 to 7 gas pressure scale 
heights. This is the first study that combines three-dimensional MHD simulations of zonal flows and 
dust particles feeling the gas pressure. The pressure bumps trap particles with St = 1 very efficiently. 
We show that St = 0.1 particles (of some centimeters in size if at 5 AU in an MMSN) reach a hundred- 
fold higher density than initially. This opens the path for particles of St = 0.1 and dust-to-gas ratio 
of 0.01 or for particles of St > 0.5 and dust-to-gas ratio 10 -4 to still reach densities that potentially 
trigger the streaming instability and thus gravoturbulent formation of planetesimals. 
Subject headings: magnetohydrodynamics (MHD) - planets and satellites: formation - protoplanetary 
disks 



1. INTRODUCTION 

Planets form as a side product in star formation. The 
general understanding on how plane t s in o ur solar sys- 
tem form was detailed in iSafronovl (jl969l ). Low-mass 
stars form out of molecular clouds which consist of 99% 
hydrogen and heli um (further re ferred to as gas), and 
1% dust and ices (|Lodder s 2003), i.e., everything that 
has a higher complexity than hydrogen molecules or he- 
lium. Those molecular clouds h ave cores between less 
than 0.1M o and more than lOM (|Krumholz et al1l2012h 
that are gravitationally unstable. Most parts of the 
mass will collapse into a newborn star. The remaining 
~ 1% of the total mass will form an acc retion disk with 
pressure supported, sub-Kepleri an gas (jWeidenschillingl 
ll977atlCassen fc Moosmanlll98l[ ) around the young star. 
Those disks have lifetimes in the order of a few million 
years (jHaisch et al.ll200H lFedele~e t al. 201 (3). Dust par - 
ticles grow due to coagulation (jWeidenschillingl Il997| ) . 
However, coagulation models show that there are sev- 
eral barriers to overcome to grow dust large enough to 
become gravitationally bound in kilometer -sized plan- 
etesimals, such as the boun cing; barrier (jZsom et al.l 
120101 : IWi ndmark et al. 2012a, fl), the fragmenta tion bar- 
rier (e.g. JBeitz et al.l l2011: Bi rnstiel et al.ll2012l, and ref- 
erences therein) , and the kilometer-size barrier (Id a et al.l 
120081 : iCuzzi e t al. 200 8|) . Dus t grow th mechanisms are 
summarized in iDominik et al.l (j2007[ ) and the review of 
IBlum fc WurmI (j2QQ8h gives an overview on the men- 
tioned barriers. 

This paper addresses the fragmentation barrier or 
meter-size barrier. Pebbles of several decimeters in size 
will drift very fast inward due to the headwind from the 
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sub-Keplerian gas (Weidenschilling 1977a). Thus, dust 
has to grow very quickly from some centimeters to sev- 
eral kilometers in size in order to avoid drifting into the 
inner region of the protoplanetary disk. 

Turbulence in protoplanetary disks around young 
stars provides pr omising mechanisms for rap id planetes- 
imal formation (IJohanse n et al.ll2007l 120111) . Shearing 
box simulations (jBrandenburg et al.l Il995[ ) are a pow- 
erful tool for analyzing the magn e torot ational instabil- 
ity (MRI; iBalbus fc Ha wlev 199lL fl998h as a source of 
turbulence. These simulations consider a local, coro- 
tating box, representi n g a sm all part of a Keplerian 
disk. iJohansen et al.l (|2009af ) reported long-lived ax- 
isymmetric sub- and super-Keplerian flows, zonal flows, 
in shearing box simulations of turbulence caused by 
the MRI. T hese zonal flows have b e en see n in several 
other local (Fromang fc Stond [2009; Stone 
20101: Simon et al. 2 Oil) and global dLvra et al.l l20j 
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Dzvurkevich et al.l 120101 : lUribe et al.l l2011t iFlock et al.l 
201 lL l2012f ) simulations using a wide variety of codes. 

Zonal flows are a product of large-scale variations 
in the magnetic field that transport momentum differ- 
entially, creating regions of slightly faster and slightly 
slower rotating gas. Large-scale pressure bumps are ex- 
cited through geostrophic balance. This creates long- 
lived over-densities that potentially trap dust particles. 
A more thorough de scription of zona l flows a nd their cre- 
ation put forward in IJohansen et al.l (j2009af ) found zonal 
flows always populating the largest radial mode avail- 
able in the local box approximation. Their largest box 
was sim ulating 10.56 press ure scale heights (H). More 
recently iSimon eFall ([2012h found a more complex struc- 
ture in their largest simulation with L x = 16 H. They 
further studied the autocorrelation function (|Guan et al.l 
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120091 ) of the magnetic field and the gas density. Both 
have a two-component structure. The first is tilted with 
respect to the azimuthal axis and highly localized. The 
second component is seen at the largest scales and can 
be associated with the (predomina n tly to roidal) back- 
ground magnetic field. iSimon et al.l (|2012f ) measure the 
radial length scale of the zonal flows to converge at 6H. 

In this paper we consider even larger physical extents 
for zonal flow structures. This gives us the opportunity 
to measure physical properties such as size and lifetime 
independent of the simulated domain. Further, we in- 
vestigate properties of the zonal flows in radially and 
azimuthally stretched boxes. We alter the radial and az- 
imuthal domain up to ~20 gas pressure scale heights. 

Add itionally, we stu dy the behavior of dust in zonal 
flows. iWhippld (fl972l ) was the first to sug gest that ax- 
i symm etric pressure bumps can trap gas. iPinilla et al.l 
(|2012f ) invoked zonal flows as a possibility to explain the 
submillimeter and millimeter-sized particles observed in 
protoplanetary disks. They used artificial static density 
bumps introduced as sinusoidal density perturbations 
with different amplitudes (e.g., A = 0.1 and A = 0.3) 
and different wavelengths (L = 0.3 . . . 3H). They found 
that a 30% density perturbation (with L = 1H) is neces- 
sary to stop the drift of the dust grains. The present work 
is the first three-dimensional MHD study that combines 
zonal flows and the reaction of dust particles on them. 

Our paper is organized as follows. In Section [2] we 
discuss the setup of the simulations in this paper. In 
Section [3] we study the zonal flow properties and their 
dependency on the physical box size. The behavior of 
dust particles in zonal flows is described in Section 2J 
A discussion and conclusions follows in Section [5] and 
Section [6] provides a summary and an outlook. 

2. SIMULATION SETUP 

We use the Pencil CodeQ a sixth-order spatial and 
third-order temporal finite difference code, for our simu- 
lations. We simulated the standard ideal MHD equations 
in a local shearing box with vertical stratification. The 
simulation boxes are centered at an arbitrary distance r 
to the star. The radial direction is denoted by x, the 
azimuthal direction by y, and the vertical direction by z. 
The Keplerian frequency is Q. We include dust particle 
dynamics, without back-reaction to the gas and without 
self-gravity. 

2.1. Gas Dynamics 

The gas velocity u relative to the Keplerian shear is 
evolved via the equation of motion 

du _ N fn\du 

~dt 



2 Qu y x — -Qu x y + Q 2 zz 

+ \j xB --VP + f u {u,p) . 



(i) 



On the left-hand side of the equation, the second and 
third terms are the advection terms by the perturbed 



velocity and by shear flow, respectively. The right-hand 
side contains the Coriolis force, the vertical component of 
the stellar gravity, the Lorentz force, the pressure gradi- 
ent, and the viscosity term. Here, Uy^ = — (S/2)Qx is the 
Keplerian orbital velocity. The magnetic field B as well 
as the current density J are calculated from the vector 
potential A using B = V x A and J = /i^Vx (V x A), 
respectively. Here, po is the vacuum pe rmeab ility. The 
viscosity term f v is explained in Section 12.2.11 

We evolve the magnetic potential with the uncurled 
induction equation 



3 A 

-dt + < >: ^ = uxB + ~2 nA y* + ^ {A) ■ (2) 

The terms on the right-hand side express the electromo- 
tive force, the stretching (creation of azimuthal magnetic 
field from radial field) by Keplerian shear and the resis- 
tivity f v (see Section E22J). 

The gas density is evolved with the continuity equation 



dA 



v dy 



-pV-u + f D (p) , (3) 



where the last term on the right-hand side describes 
mass diffusion (see Section I2.2.3|) . We use an isother- 
mal equation of state P = c 2 p, where the speed of sound 
is c s = HQ; H is the gas pressure scale height. 

2.2. Dissipation 

Maxwell and Reynolds stresses as well as the MRI re- 
lease kinetic and magnetic energy at large scales. This 
energy cascades down to small scales. Since numerical 
simulations have a finite resolution, this small-scale en- 
ergy needs to be dissipated. We use numerical diss ipatio n 
in the form of hyper- and shock viscosity (Section T2. 2. 1|) . 
hyper-resistivity (Section I2.2.2p . and hyper- and shock 



diffusion (Section 12.2. 

2.2.1. Viscosity 
The viscosity term f v in Equation (pQ) is expressed by 



(S (3) ■ Vlnp 

+ ^sh [ VV • u + (V • u) (V • In p)\ 

+ (V^ sh )V-u. (4) 

We restricted our models to hyper- (^3) and shock (v s h) 
viscosity. Thus, the regular Navier-Stokes viscosity term 
is neglected. The third-order rate-of-strain tensor is 
defined by 

,(3) _ 9 b Uj 



dx b 



(5) 



The high-order Laplacian V 6 in Equation (J2J) is ex- 
panded as V 6 = d 6 /dx 6 + d 6 /dy 6 + d 6 /dz 6 . Further- 
more, the shock viscosity is expressed by 

^sh = c s h (max [—V • u] + ) min (5x, 5y, Sz) 2 . (6) 

In the fashion of I von Neumann fc Richtmverl (pj950) it is 
proportional to positiv43 flow convergence. We take the 



1 Details on the Pencil Code and download information can be 
found at http://www.nordita.org/software/pencil-code/ 



2 Symbolized by the plus sign in Equation (|6j). We only apply 
shock viscosity where the velocity flow is converging. 



Planetesimal Formation in Zonal Flows 



3 



maximum over five zones, and smoothed it to the sec- 
ond o rder. As suggested bv Ivon Neumann fc Richtmverl 
(jl95Qh . we set the shock viscosity coefficient to c s h = 1.0 
to dissipate energy in shocks at high z above the mid- 
plane of the disk. 

2.2.2. Resitivity 

The effects of resistivity are captured by the term 

f v = V3 V 6 A, (7) 
where 773 is the hyper-resistivity. 

2.2.3. Diffusion 
Mass diffusion is computed with 

f D = D 3 V 6 p + D sh V 2 p + VAh • Vp , (8) 

where D3 is the hyper-diffusion parameter and D s h is 
expanded as in Equation (j6]). 

2.3. Dust Dynamics 

Dust particles are simulated as individual super- 
particles i with position X{ and velocity v\. Each super- 
particle position is evolved with 

dx^ 

(9) 



dt 



The change of velocity for each particle is evolved 
through 



dt 



--2ttv^x- 



1 

Tf 



1 

~ 2 
u(x^) 



y - Vtrzz 



(10) 



where the first and second terms are due to the Coriolis 
force. The third term corresponds to the vertical gravity 
of the star. Particles only feel the gas drag (the last term 
in Equation ([TO]) ) of nearby cells, but are not subjected 
to pressure or Lorentz forces. Tf denotes the friction time, 
a measure for the size of the particles. 

2.4. Boundary Conditions 

For our simulations, we use shearing box boundary 
conditions in radial (shear-periodic) and azimuthal (pe- 
riodic) directions. In the vertical direction we also use 
periodic boundary conditions. Although periodic bound- 
ary conditions in vertical direction are not physical, these 
boundary conditions conserve the average flux of the 
magnetic field. Simulations with outflow boundaries (not 
included in this paper) showed no considerable mass flux 
across the vertical boundary and did not change the av- 
erage properties of the zonal flow. 

2.5. Dimensions 

We use the dimensionless unit system c s = Q = po = 
po = 1. Velocity is measured in units of the local sound 
speed c s . Gas velocities are always denoted by u whereas 
particle velocities are always denoted by v. All velocities 
are differences to the Keplerian orbital velocity vk = 

(0,i4°\o), where u^P = —(3/2)Qx. Time is measured 
in units of the local orbital time T or b = 2irQ~ 1 . Length 



measures are in units of the pressure scale height H = 
c s £7 -1 . Density is stated in units of the initial mid-plane 
gas density po. Magnetic field strength is measured in 
units of c s (/ioPo) _1 - Energy and stress are in units of the 
mean thermal pressure in the box (P) = c 2 (p). 

Since our simulations are dimensionless, they can be 
placed at any distance r to the star. Only by defining a 
global pressure gradient dP g i b a i/dr, which balances the 
Coriolis force in 



1 <9Pgi bal 

P 



dr 



29, Av, 



(11) 



we restrict our simulations to a specific distance to the 
star where the chosen pressure gradient applies. The pa- 
rameter Av = v!y^ — Uy is the difference to the azimuthal 
Keplerian velocity. We fix Av = 0.05c s (see also Sec- 
tion 12. 6p . Numerically, the global pressure gradient acts 
as an external force on gas and dust. 

2.6. Initial Conditions 

The gas density is set to an isothermal hydrostatic 
equilibrium p{z) = po exp (— z 2 /2H 2 ). We start with 
random noise fluctuations in the gas velocity with 



Su = 10- 



Cs. 



The azimuthal component of the 



magnetic vector potential is initialized with A y = 

Aq cos (k x x) cos (k y y) cos (k z z) where throughout k x = 



k z 



4.76H- 1 and A = 0.04c s (/i po) _1 . 



Particles are released after the gas turbulence is satu- 
rated. We measured this to be after 20T or b for the largest 
runs. For convenience, we used the same saturation time 
for all our simulations. Particles have a Stokes number 
of St = Tffi = 1, unless otherwise stated. The initial par- 
ticle distribution is Gaussian in z and uniform in x and 
y. The parti cle velocity is initialize d with the station- 
ary solution (Nakaga wa et al.l |1986[ ) for the radial and 
azimuthal velocity 



Vx_ 

Cs 

Vy 



2Av 



Tfft + (Tfft)- 1 

Av 



1 + (Tff2) 2 " 

We get Av from the solution of Equation JTT} 



Av 

C s 



■-(- 

2\r 



d\nP 
9 In r 



(12) 



(13) 



We initialized Av = 0.05c s for our simulations. 

2.7. Simulation Parameters 

The parameter space covered by our simulations 
is summarized in Figure [U The vertical extent 
is always set to L z = 2.64i^Q One simula- 
tion set (A) covers the boxes with a squared base, 
i.e., radial and azimuthal extent are kept the same: 
L x = L y = {1.32, 2.64, 5.28, 10.56, 21.12}#. 
These are marked with blue boxes in Figure [1] and 
are called runs 5, M, L, XL, and XXL. The devia- 
tion to the global density profile in the largest box 

3 L = 1.32 has been chosen as the basic box size, because 
L x = 1.32 approximately mark s the transition from s ubsonic to 
supersonic Keplerian shear flow ( Johansen et al. 2009a). 
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Table 1 

Run Parameters 



Simulation Set 


Run 


L x x Ly x L z 


N x x N y x N z 


^3 




m - 


~-D 3 


^particles 


St 


Shear 


At 


(1) 


(2) 


(3) 


(4) 






(5) 




(6) 


(7) 


(8) 


(9) 


A 


S 


1.32 x 1.32 x 2.64 


36 x 36 x 72 


4.0 


X 


10" 


10 


62,500 


1.0 


FDA 


121 


A,B,C 


M 


2.64 x 2.64 x 2.64 


72 x 72 x 72 


4.0 


X 


io- 


10 


250,000 


1.0 


FDA 


121 


A,E 


L 


5.28 x 5.28 x 2.64 


144 x 144 x 72 


4.0 


X 


io- 


10 


1,000,000 


1.0 


FDA 


121 


A, hi 


XL 


ID. 5b X lU.oo x 2.64 


288 X 288 X 72 


4.0 


X 


10~ 


10 


4,000,000 


1.0 


FDA 


121 


A 


XXL 


21A2 X 21.12 X 2.64 


57b x 57b x 72 


4.0 


X 


10~ 


10 


A nnr\ nr\n 

4,000,000 


1.0 


FDA 


121 


B 


x-S 


1.32 x 2.64 x 2.64 


36 x 72 x 72 


4.0 


X 


10" 


10 


125,000 


1.0 


FDA 


121 


B 


x-L 


5.28 x 2.64 x 2.64 


144 x 72 x 72 


4.0 


X 


10" 




500,000 


1.0 


FDA 


121 


B 


x-XL 


10.56 X 2.64 X 2.64 


OOO v , TO v , TO 

288 X 72 X 72 


4.0 


X 


10~ 


10 


1,000,000 


1.0 


FDA 


121 


C 


y-s 


2.64 x 1.32 x 2.64 


72 x 36 x 72 


4.0 


X 


io- 


10 


125,000 


1.0 


FDA 


121 


C 


y-L 


2.64 x 5.28 x 2.64 


72 x 144 x 72 


4.0 


X 


10" 


10 


500,000 


1.0 


FDA 


121 


C 


y-XL 


2.64 x 10.56 x 2.64 


72 x 288 x 72 


4.0 


X 


10" 


10 


1,000,000 


1.0 


FDA 


121 


D 


LspecMR 


5.28 x 5.28 x 2.64 


144 x 144 x 72 


4.0 


X 


10" 


10 


1,200,000 


0.01... 100 


FDA 


121 


D 


LspecHR 


5.28 x 5.28 x 2.64 


256 x 256 x 128 


2.0 


X 


10" 


11 


120,000,000 


0.01... 100 


FDA 


121 


D 


LspecMRs 


5.28 x 5.28 x 2.64 


144 x 144 x 72 


4.0 


X 


10" 


10 


14,000,000 


0.01... 1.0 


FDA 


121 


D 


MspecMRb 


2.64 x 2.64 x 2.64 


72 x 72 x 72 


4.0 


X 


10" 


10 


6,000,000 


1.0... 100 


FDA 


223 


E 


L_SAFI 


5.28 x 5.28 x 2.64 


144 x 144 x 72 


4.0 


X 


10" 


10 


1,000,000 


1.0 


SAFI 


121 


E 


XL.SAFI 


10.56 x 10.56 x 2.64 


288 x 288 x 72 


4.0 


X 


10" 


10 


4,000,000 


1.0 


SAFI 


121 



Notes. Column 1: simulation 
resolution. Column 5: dissipation 



set. Column 2: name of run. Column 3: box size in units of pressure scale heights. Column 4: grid 
coefficients. Column 6: number of particles in simulation. Column 7: Stokes number St = Tff2. Column 
8: shear advection scheme. Column 9: total run time in orbits T or b. 



21.12 



10.56 



5.28 



2.64 



1.32 











XXL 

576x576x72 
Af p =4xl0 6 




y-XL 

72x288x72 
N=10 6 




XL 

288x288x72 
Af=4xl0 6 






y-L 

72x144x72 
Af p =5xl0 5 


L 

144x144x72 
N=10 6 




x-S 

36x72x72 
N=l 25,000 


M 

72x72x72 
N p =250,000 


x-L 

144x72x72 
AT p =5xl0 5 


x-XL 

288x72x72 
Af p =10 6 




S 

36x36x72 
Af=62,500 


y-S 

72x36x72 
# p =125,000 







1.32 2.64 



5.28 
L/H 



10.56 21.12 



Figure 1. Parameter space of radial and azimuthal box sizes that 
was simulated for this paper. Every simulation has a vertical extent 
of 2.64H . The first line in each box states the name of the run, 
the second line the number of grid cells used, and the third line 
gives the number of simulated super-particles. More details on all 
simulations are found in Tabled 



can be quite severe at the inner and outer bound- 
ary of the largest simulation. Thus, the results from 
run XXL have to be treated with caution. Another 
set of simulations (B) varies the radial size of the box, 
L x = {1.32, 2.64, 5.28, 10.56} H, with constant box 
size in azimuthal direction, L y = 2.64H. This set is 
marked red in Figure [1] and includes runs x-S, M, x-L, 
and x-XL. The third set of simulations (C) varies the 
azimuthal extent, L y = {1.32, 2.64, 5.28, 10.56} if, 
while the radial extent is kept constant, L x = 2.64H. 
This set includes runs y-S, M, y-L, and y-XL (marked 



yellow in Figure [T]). All simulations are stratified and 
have dust particles with different couplings to the gas. 
The simulations displayed in Figure [1] have particles with 
a Stokes number of St = 1. 

Details on run parameters of those and six more sim- 
ulations are found in Table [TJ The first set of simula- 
tions (A,B, and C) in Table [1] are the simulations with 
medium resolution, i.e., 36 grid cellfl per 1.32 pressure 
scale heights. Simulation set D was carried out to in- 
vestigate the behavior of different particle sizes in the 
presence of zonal flows. Run LspecMR is very much like 
run L, but with 12 different particle Stokes numbers. The 
run LspecHR has a resolution of 64 grid cells per 1.32i7. 
Runs LspecMR and LspecHR have 12 different parti- 
cle species, with Stokes numbers of St = 0.01 . . . 100.0. 
The runs LspecMRs and MspecMRb have particles with 
Stokes numbers of St = 0.01 ... 1.0 and St = 1.0 .. . 100.0, 
respectively. These two simulations were carried out to 
study particle behavior with more particles per grid cel0 
at medium resolution. The corresponding sizes for differ- 
ent protoplanetary disk models are found in Section 15.21 

Simulation set E is a comparison of runs L and 
XL to the same runs (L_SAFI and XL.SAFI) with 
the Shear Advection by Fourier Interpolation (SAFI) 
scheme. Here, all variables q(x,y, z) are transformed into 
Fourier space in the ^-direction to get q(x,k y ,z). Then 

each Fourier mode is multiplied by exp [ik y u y \x)St\ to 

shift by Uy\x)8t in real space and is inverse Fourier 
transformed to real space. This method reduces the ad- 
vection error to the standard Finite Difference Advection 
(FDA) scheme in the P encil Code (more deta ils on FDA 
and SAFI are found in I Johansen et al"]l2009a| ). 

4 We chose 36 grid cells instead of the usual 32 grid cells. That 
choice was done due to the architecture (12 CPUs per node) of the 
used cluster, THEO in the MPG computing center in Garching. 

5 ~ 9 and ~ 16 particles per grid cell for runs LspecMRs 
and MspecMRb compared to ~ 0.8 particles per grid cell for 
run LspecMR. 
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Table 2 

Turbulence Properties 



Run 


( 


H) 




( 






( 


H) 






( 








( 






( 






(pu x Uy 


> 


(- 


B X By) 


a 




(1) 




(2) 






(3) 






(4) 








(5) 








(6) 






(7) 






(8) 






(9) 


(10) 




S 


2.1 


X 


10" 


-3 


3.3 


X 


10" 


-3 


1.5 


X 


10" 


-3 


8 


.8 


X 


10" 


-4 


6, 


.4 


X 


10" 


-3 


3.5 


X 


10" 


-4 


7.7 


X 


10 


-4 


3.4 


X 


10~ 3 


2.8 x 10- 


-3 


M 


3.9 


X 


10" 


-3 


5.2 


X 


10" 


-3 


2.2 


X 


10" 


-3 


1 


.9 


X 


10" 


-3 


1 


.2 


X 


10" 


-2 


7.6 


X 


10" 


-4 


1.6 


X 


10 


-3 


6.6 


X 


10- 3 


5.5 x 10" 


-3 


L 


5.0 


X 


10" 


-3 


5.6 


X 


10" 


-3 


2.3 


X 


10" 


-3 


2 


.0 


X 


10" 


-3 


1 


.3 


X 


10" 


-2 


8.0 


X 


10" 


-4 


1.9 


X 


10 


-3 


6.9 


X 


10- 3 


5.9 x 10" 


-3 


XL 


5.2 


X 


10" 


-3 


5.0 


X 


10" 


-3 


2.1 


X 


10" 


-3 


1 


.7 


X 


10" 


-3 


1 
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run. Columns 2-4: kinetic e nerg y. Columns 5-7: magnetic energy. Column 8: 
on- value, following Equation (|14|) . Stresses and energies have been normalized 
in the box, (P) = c*{p). 



Notes. Column 1: name of 
Maxwell stress. Column 10: i 



Reynolds stress. Column 9: 
to the mean thermal pressure 



Every simulation is run for 121 local orbits T or b = 
27rf2 -1 , except for LspecMRb which runs for 223T or b in 
order to follow the evolution of the slowly settling large 
particles. After 20T or b, when the initial conditions are 
sufficiently forgotten and the turbulence saturated, the 
particles are started. 



3. ZONAL FLOW PROPERTIES 

Turbulence properties are summarized in Table O The 
kinetic and magnetic energy as well as the Reynolds and 
Maxwell stress almost doubles when increasing the box 
size from (1.32#) 2 x 2MH (run S) to (2.64#) 3 (run M). 
Further increasing the box size does not change the re- 
sulting energies and stresses by much. The radially short 
box of run x-S with 1.32H x (2.64H) 2 has similar results 
on these values. However, the azimuthally short box of 
run y-S has turbulent energies and stresses comparable 
to run S. These measurements show that the turbulence 
parameters are saturated for boxes with an azimuthal 
extent of at least 2 .64H . This confirms the results from 
Fr omang fc Stone (2009) who found that the turbulence 
properties do not change when the box size is increased 
radially, if the azimuthal dimension is large enough. The 
a- value ( Sh akura fc Sunvaevlll973l ) in Column 10 in Ta- 
ble [2] is calculated via 



iSimon et all (|2012f ): 



(15) 



2 ((pU X Uy) 

3 



(B X By)) 



(P) 



(14) 



where (P) = c 2 s (p). The factor of 2/3 originates from 
the shear parameter q = — dhiQ/dhiR. We use q = 3/2, 
appropriate for a Kepleri an disk. For further details see 
iBrandenburg et all (|1995l page 748). The Maxwell stress 
is around three times higher than the Reynolds stress and 
thus dominates the a- value. 

In order to verify that our numerical resolution is suf- 
ficient, we examined the quality factor as described in 



where the Alfven speed is defined as |^ a ,j| 2 = (Bj) 2 /(p)- 
The notation ( x) denotes volume av eraging, x shows a 
time average. iSorathia et al.l (|2012f ) show that Q z > 
10 — 15 for poorly resolved azimuthal quality factors 
(Q y ~ 10) is required to resolve the MRI. Larger val- 
ues of the azimuthal quality factor (Q y > 25) allow for 
lower vertical quality factors. The azimuthal component 
of the magnetic field is very well resolved (Q y > 25) for 
all simulations, but runs S and y-S. The vertical com- 
ponent has values between 6 and 8. We thus conclude 
that all simulations, but runs S and y-S have sufficient 
resolution for the MRI. 

In Figure [2] a snapshot of the runs y-XL, XL, M, and 
x-XL are shown in scale, giving a real size comparison 
of high- and low-pressure regions. The large-scale sinu- 
soidal form of the dominant mode is observable in these 
plots. The higher modes are much shorter lived and seem 
to be non-axisymmetric density wave s affected by the 
shear ([Heinemann fc Papaloizoul 1200 9). The amplitudes 
of the pressure differences are higher in azimuthally large 
boxes. Only the axisymmetric density waves are long- 
lived and strong enough to make up a significant contri- 
bution to the pressure bump structure in an azimuthal 
as well as a temporal average over some local orbits. The 
azimuthal average of Figure [2] is seen in Figure [3l Here, 
the axisymmetric structure is clearly visible and a strong 
correlation between the particle location and a positive 
radial gradient of the gas density is seen. The black dots 
in Figures [2] and [3] show the radial and azimuthal posi- 
tion of every 100th particle. The particles are trapped 
by the axisymmetric pr essure bumps. Also, the shapes 
of spiral density waves (jHeinemann fc Papaloizoul [2009) 
can be seen in the structures. The particle distribution 
with respect to the gas flow will be discussed in more 
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Figure 2. Collage of four gas surface density representations of the runs y-XL, XL, M, and x-XL. Each snapshot was taken after 85T or £,. 
These plots show that gas overdensities are most pronounced in the largest box. The non-axisymmetric structures have very short lifetimes; 
less than a tenth of an orbit. The pressure bump structures are more visible when the density is averaged over the azimuthal direction 
(Figure [3]). The black dots represent the position of every 100th particle, integr ated in vertical direct i on. T he particles are trapped both 
in axisymmetric pressure bumps and in spiral density waves as described in Heinemann &; Papaloizou (2009). 



detail in Section [H 

All of our simulations show signs of zonal flows. 
Strength and lifetime of the zonal flows and the asso- 
ciated pressure bumps differ very much with the physi- 
cal box size. Space-time plots of all different simulation 
sizes are shown in Figure [H and the upper left panel of 
Figure [5l The pressure bumps are generally more pro- 
nounced in simulations with a larger radial extent. Sim- 
ulation set A strictly follows this general trend. Pressure 
bump features grow in strength and lifetime with the 
physical box size, always staying at the largest radial 
scale. This rule applies to all but the largest runs XL, 
x-XL, and XXL. There, instead of the formerly predom- 
inant k x = 1 (uu x = 2i\k x jL x in Fourier mode s'm (uj x x)) 
mode, the mode k x = 2 (higher modes for run XXL) is 
occupied by the pressure bumps. In simulation set B, 
the strength and size of the pressure bumps converges 
for simulations with a radial extent larger than 5.2SH. 
The lifetime of the pressure bumps even decreases for the 
largest simulation in this set. Simulation set C is qualita- 
tively different from the other simulation sets. Strength, 
size, and lifetime of the pressure bumps seem to be in- 
versely proportional to the azimuthal extent of the box, 
when the vertical and radial box size s are kept constant . 
Thi s effect was a lready seen in, e.g.. iSimon eFall (|2Q12h 

d lFlock et all (|2012f ). Both groups show that the mag- 
netic field consists of two components: a local turbulent 
component that is responsible for the zonal flows and 
a global azimuthal component. Since the total energy 
stays approximately constant, the local component gets 
weaker and consequently zonal flows as well as axisym- 



metric pressure bumps get weaker too. 

Figure [5] additionally shows space-time plots of the az- 
imuthal gas velocity and the radial gradients of gas den- 
sity and azimuthal gas velocity of run XXL. In the right 
panels, the position of the highest dust density are shown 
as dots. Particles get clearly slowed by the maxima of 
the azimuthal gas velocity, i.e., the large-scale maxima 
of the pressure gradients. The velocity has large-scale 
structures that are very similar to those of the density 
gradient, as expected in geostrophic balance. Thus, the 
structure of the velocity gradient can be approximated as 
the large-scale structure of the second derivative of the 
gas density. It is shown in the lower right panel that the 
particles get stopped at the minima of the radial deriva- 
tive of the azimuthal gas velocity (and thus at minima of 
the second derivati ve of the gas densit y) as analytically 
predicted (see, e.g. JKlahr fe Linl 120011 ). 

We calculated the correlation time of the pressure 
bumps and t he zonal flows in the s ame way as it was 
calculated in Uohansen et al.l (j2009aj ) . We use the den- 
sity p, averaged over azimuthal and vertical directions, at 
a given time t. Then, we average over each point in radial 
direction the time it takes for the density at each point 
to change by a value corresponding to the standard devi- 
ation of the gas density. These measurements are taken 
for every local orbit. The measurements are averaged 
over the time between saturation of the turbulence and 
a time when the correlation does not extend the corre- 
lation time to the final time of the simulation. Finally, 
the averages are multiplied by two, in order to cover the 
full temporal extent of the correlated structures. The 
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x/H xlH 

Figure 3. Surface density distribution of Figure [2] averaged in azimuthal direction and averaged over the mean surface density. This 
reveals axisymmetric pressure bumps and valleys. Particles are trapped on the inner side of the density maxima, at places with a positive 
density gradient to overcome the negative global pressure gradient. These pressure bumps are stable for many orbits (compare Figures [4] 
and [5}. 



correlation times measured in this fashion are in good 
agreement with the lifetime of the overdensities that is 
seen in Figures |4] and [5] However, a change of position of 
the structures, as seen in run XL (check Figure 2]), is not 
accounted for. Thus, correlation times are more likely to 
be underestimated than overestimated. Also, we cannot 
be entirely sure whether this behavior is really drift or 
structure decay and reformation. 

The results of the correlation time determination are 
shown in Table [3] and in the upper panel of Figure [6j For 
the diagonal simulation set, (A), the correlation time in- 
creases with box size. It seems to saturate toward the 
largest box size. The trend to longer correlation times is 
also evident for simulation set B. Here only run x-XL has 
a shorter correlation time than expected. This might be 
an effect of the strongly stretched simulation box. The 
correlation time decreases slightly with an increasing az- 
imuthal box size in simulation set C (not shown in the 
figure). The lower panel in Figure [6] shows a measure- 
ment for the physical size of the zonal flow features. We 
Fourier-transformed the vertically and azimuthally av- 
eraged gas density and azimuthal gas velocity for each 
time step and averaged the amplitudes of the first four 
modes over the time of 20 . . . 120 local orbits. The length 
was normalized for the size of the simulation box, to get 
the physical size of the modes by \ x = L x /k x with the 
wave number k x . The turbulence is always strongest at 
the largest modes for simulations with L x < 5H. The 
highest amplitude for both quantities in the largest sim- 
ulation domain is found between 5H and 7H (up to 10H 
for p). These measurements are also found in Table [3l 



The runs L_SAFI and XLSAFI were carried out to 
compare the turbulence and zonal flow parameters with 
the runs L and XL. They were run to check that zonal 
flows are no effects from the shear advection scheme that 
was used in the Pencil Code. Comparing the values in 
Tables [2] and [3] shows that there is little change in the 
measured properties of the zonal flows and the associ- 
ated pressure bumps. However, the computation time 
increases if one uses the SAFI scheme. Thus, this scheme 
was only used to confirm our results. 

4. PARTICLE BEHAVIOR IN ZONAL FLOWS 

Particle accumulations and planetesimal formation can 
occur in clumps and filaments of the overdensities in the 
dust. In our simulations, we do not include gravitational 
interaction between the particles. Thus, we only study 
the passively developed overdensities of the dust to see 
when and whether overdensities sufficient for the stream- 
ing instability can be reached. By not having explicit 
feedback one can retroactively study the concentration 
factor for various initial dust-to-gas ratios. Simulations 
including feedback will have to be done in future stud- 
ies. Figure [2] shows the position of every 100th parti- 
cle in selected simulations. These plots clearly show the 
trend for particles to accumulate in the downstream of 
high-pressure regions. Particles are pulle d toward pres- 
sure gradient maxima (jKlahr fc Linll200Th . In the upper 
right panel in Figure [2j a snapshot of run XL after 85 
local orbits is shown. The particles clump up at posi- 
tions just left of the maxima in the k x = 1 of the gas 
density; these are the locations of positive zonal flows, 
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Zonal Flow Properties 
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Notes. Column 1: name of run. Column 2: root-mean-square density p rms = y ((p — ~p) 2 )- Columns 3-5: Fourier amplitude of radial 
density modes k x = 1 ... 3, normalized by mean density in the box. Columns 6-8: Fourier amplitude of azimuthal velocity modes 
k x = 1 ... 3 with u y = u y — u y . Column 9: correlation time, in orbits T = 27rf2 -1 , of the largest radial density mode. 



i.e., regions where the azimuthal gas velocity is higher 
than the pressure-supported Keplerian flow. 

4.1. Particles in Zonal Flows 

In the upper right panel of Figure the azimuthal gas 
velocity development of run XXL is shown, overplotted 
with the position of the most massive clump for each 
time step. The azimuthal gas velocity coincides with 
the derivative of the gas density, but it is much eas- 
ier to interpret. The speckled structure of the deriva- 
tives comes due to the high power in the smaller scales. 
However, the large-scale structure is still visible and the 
geostrophic correlation between the structures of u y (x, t) 
and d/dx[p(x,t)\ is directly observed. Since they have 
the same large-scale structure, the particle position is 
much easier interpreted at the azimuthal gas velocity plot 
than on the density gradient plot. Sometimes the radial 
displacement from one orbit to the next is too large to 
be explained by radial drift. That happens when an- 
other clump becomes more massive than the previous 
one. These particles accumulate in regions with high az- 
imuthal gas velocities (see upper right panel in Figure [5]). 
The only time when this is not true is at times from 80 
to 100 local orbits. In this period, an inward-drifting 
clump stayed coherent during the time of its drift. The 
drift velocity of the most massive particle clump is indi- 
rectly encrypted in this plot. Particles are drifting much 
slower when they are trapped by a pressure gradient. As 
all particles drift inward this leads to accumulation of 
particles in regions where the perturbed pressure gradi- 
ent is positive. 

The maximal accumulation of particles for runs XL 
and y-S are plotted in the top panel in Figure 
The second panel shows the evolution of the quantity 

(u y — (u y )y Z )l, a measure for the strength of the zonal 
flows. The third panel in Figure [71 shows the evolu- 
tion of the strength of the gas density enhancement as 
yj (p — (p) yz )x- Comparing the second and third panels, 



one can see a clear correlation between the zonal flow 
strength and the gas density enhancements. The bottom 
panel in Figure [71 shows the evolution of the a-parameter, 
calculated as in Equation ([T4]) . 

The maximum of the dust overdensity that occurs dur- 
ing one simulation is plotted against the box size in the 
upper panel of Figure [U The general trend shows that 
radially larger boxes have higher particle concentrations. 
An increased azimuthal extent does not have an effect on 
the particle concentrations. The most surprising result 
is in run y-S. It shows a very high particle concentration 
that occurs early in the simulation (compare Figure [7j). 
This is most likely a stochastic coincidence. The lower 
panel in Figure [8] shows a plot of the maximum dust over- 
density against the correlation time of the zonal flows. 
The error margin are calculated with the standard de- 
viation of the temporal evolution of the two quantities. 
We see a clear trend that denser particle accumulations 
develop with longer correlation times. The distribution 
can be fitted by a power law. This gives an exponent of 
dlogp/<ilogT C orr = 0.38 ± 0.05. The one point that does 
not overlap with the error margins of the fit is from run y- 
S. If we take the maximum of the top panel in Figure [7] 
after the two first maxima (i.e., after 45T or b) and plot 
this value again in the parameter space of Figure [8j we 
get the position marked with the blue square. It agrees 
well with the error margins of the fit. 

In isothermal geostrophic balance, 2pQu y = c^dp/dx, 
the azimuthal gas velocity follows the radial density gra- 
dient. That this is true for large scales as shown in Fig- 
ure [U The upper left panel shows the evolution of the 
azimuthally and vertically averaged azimuthal compo- 
nent of the gas velocity. Overplotted are the locations of 
the maxima in the dust density. In the upper right panel 
the dust density evolution of the same run L is plotted. 
In comparing the location and times of the maxima and 
minima on these two plots, one clearly sees that max- 
ima in the dust density occur often at times and loca- 
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Figure 4. Evolution of the gas density perturbation of all runs from simulation sets A, B, and C. Run XXL is shown in Figure GD The 
density is averaged in vertical and azimuthal direction and plotted in radial direction over time. The lifetime, the size, as well as the 
strength of the pressure bumps are clearly increasing with increasing box size in simulation set A, i.e., runs S, M, L, XL, and XXL. In 
simulation set B, i.e., runs x-S, M, x-L, and x-XL, we have the same increase of lifetime, size, and strength of the pressure bumps. Only 
for the very large simulation (L x = 10.56H), there is no apparent difference in pressure bump size and strength to run x-L. For simulation 
set C, i.e., runs y-S, M, y-L, and y-XL, the strength of the pressure bumps is apparently constant throughout this set of simulations. Even 
the lifetime decreases slightly with increasing box size. 



tions where one finds maxima in the gas velocity. Two 
attempts to quantify this observation are shown in the 
lower row of Figure [9l In the left panel, the particle 
density and the azimuthal velocity from the two upper 
panels are plotted against each other, regardless of po- 
sition and time. In the right panel, a snapshot of the 
simulation (as in Figure [2]) was taken at 85 local orbits, 
the time when the maximum dust density enhancement 
occurs. The particle density as well as the azimuthal gas 
velocity were integrated in vertical direction and plotted 
against each other, regardless of their radial or azimuthal 
position in the simulation, in this scatter plot. In order 
to visualize high densities of points in these plots, we 
computed a two-dimensional histogram of the scattered 



points. This is indicated by the color scale, showing the 
amount of points in each of the boxes in the scatter plot 
space. There is a clear trend for high dust density con- 
centrations to appear at high gas velocities. Without 
radial drift particles would concentrate where u y = 0, 
i.e., between the sub- and super-Keplerian flow. Due 
to the radial drift particles accumulate slightly down- 
stream at the formed pressure bumps. Those happen to 
be at the maxima of the azimuthal gas velocity. With the 
geostrophic balance, high velocities are also regions of a 
high radial density gradient. These plots prove that the 
particles in the simulations are trapped by the long-lived 
pressure gradients that occur due to stable zonal flows. 
If the dust-to-gas ratio increases to values larger than 
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Figure 5. Top panels show the evolution of the gas density perturbation and the azimuthal gas velocity of run XXL, whereas the bottom 
row shows the radial derivative of these quantities. The derivatives are very speckled, since small-scale fluctuations give stronger amplitudes 
to the derivatives. However, the underlying large-scale structure is still visible. The azimuthal gas velocity follows the radial gas density 
gradient, as expected for a geostrophic balance. Hence, it is possible to interpret the radial derivative of the azimuthal gas velocity as the 
second derivative of the gas density. In the upper right panel, the black dots represent the position of the most massive particle clump in 
the simulation at each time. It is clearly shown that particles get trapped in regions of positive zonal flow downstream of pressure bumps. 



unity , the streaming inst a bility (lYoudin fc Goodman 
2005; Uohansen fc Youdinl 120071 : lYoudin fc Johansen 
20071 ) is triggered. This increases the dust density further 
on timescales shorter than an orbital period. To follow 
the streaming instability development, the back-reaction 
of the dust particles to the gas phase must be considered 
in future numerical simulation. This effect was neglected 
in this set of simulations. Otherwise the initial dust-to- 
gas ratio would have been an additional free parameter 
to be studied. 

4.2. Radial Drift 

Radial drift velocities of the particles in the simula- 
tions with different box sizes are shown in Figure [101 
The upper panel shows the measured and expected ra- 
dial drift of two simulations (M and XL). They show 
that particles drift slower in turbulent simulations than 



they would in a laminar disk. However, the size of the 
simulation has little effect on the actual drift velocity, 
as shown in the lower panel of Figure [TUJ It shows a 
time average of the particle drift velocity plotted against 
the box size. The uncertainties are too large to reveal 
a trend. Thus, the reduction of the radial drift velocity 
apparently only depends on the amplitude of the zonal 
flow, but not on the correlation time. Looking at the 
largest run XXL, we can estimate that the radial drift 
gets reduced by about 28% (drop of the absolute value 
from 0.05c s to (0.036 ± 0.003)c s ). 



4.3. Clustering 

The clustering degree of the particle distribution can 
be estimate d with the distr ibution of the dust surface 
density S p (|Pan et al .11201 If ). The initial distribution is 
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Figure 6. Upper plot shows the correlation times of the largest 
radial density mode against the radial box size. The lines corre- 
spond to the simulation sets A and B. The results from simulation 
set C are omitted for visibility. The correlation time r CO rr grows 
for boxes with a larger radial extent. Only run x-XL does not fol- 
low this trend. The large ratio L x /L y may prohibit formation of 
stable zonal flows. The two lower plots show the first four am- 
plitudes of the radial Fourier modes of the gas density and the 
azimuthal gas velocity against their real size \x — Lx I kx 5 kx is 
the wave number of the corresponding Fourier mode, defined by 
uj x = ^Trkx/Lx in the Fourier mode sin (oj x x). The lines connect 
the amplitudes of different Fourier modes for one simulation. Both 
quantities have most of their power in the largest modes. Only in 
the largest simulations, the power in the largest modes decreases. 
There the maximum is between 5H and 7H. 

represented by a Poisson distribution (see Figure [TT])FI 
For this plot, we binned the measured dust surface den- 
sity of a snapshot. We then normalized them to the 
amount of grid cells. About three local orbits after the 
particles feel the gas drag, the shape of the distribution 
function is saturated. We averaged the distribution over 
the time of 23 . . . 121T or b. We see at the high density end 
of the distribution that higher densities develop in larger 
boxes due to the higher number of available particles. 
Thus, the clustering properties do not depend strongly 
on the strength or lifetime of the zonal flows (compare 
Figure [8j bottom). 

4.4. Different Particle Sizes 

So far we only considered simulations with one parti- 
cle species, i.e., St = TfQ = 1. We take the simulation 

6 Run XXL was not included in this figure, because the number 
of particles per grid cell was different to the other runs. 




0.001 1 1 

20 40 60 80 100 120 

Figure 7. Time series of runs XL and y-S. The plots show (from 
top to bottom) the maximum of the dust density, the root mean 
square of the azimuthal gas velocity and the gas density, and the 
a- value (Equation (|14[0 . The dust overdensities of run XL have 
a higher base than those of run y-S. The latter has some spikes 
in the beginning, but is lower for most time of the simulation. 
The two panels in the middle show that the azimuthal gas velocity 
and the gas density are correlated. Both plots show maxima and 
minima at the same time, while a is rather stable with time. The 
time-averaged a- values for all simulations can be found in Table [2] 

size that simulates one fully extended zonal flow and in- 
vestigate 12 different particle species. The particle sizes 
range from St = 0.01 to St = 100. We choose run L with 
the dimensions 5.28H x 5.28H x 2.64H as simulation size 
for the last simulation set. For one simulation we used 
a smaller box, because the integration time had to be 
increased be a factor of two to give the particles with 
the high Stokes numbers the opportunity to react on the 
pressure differences. 

4.4.1. Drift Velocity and Particle Densities 

The results are shown in Figure [T21 The upper left 
panel shows the negative of the radial velocity of the par- 
ticles, averaged over all particles of a certain size and over 
time. The four different simulations match very well. 
The plot shows that particles with St = 1 drift fastest 
inward, also with turbulence in the simulations. On 
both sides the inward drift velocity decreases with simi- 
lar slopes. The key to the different colors and symbols is 
in the lower right panel. Overplotted, in a dashed gray 
line, we find the analytical prediction (following Equation 
(|T2|) ) for the radial drift in a laminar disk. The difference 
to the prediction is shown in the lower sub-panel. Large 
particles generally drift slower according to the steady- 
state solution and their coupling to the gas is also much 
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Figure 8. In the upper panel, the highest peak in the time series 
of the dust overdensity (the top panel in Figure[7|) is plotted against 
the size of the simulation box. The diagonal simulation set A is 
marked by the blue line. The dust overdensity increases with box 
size, until it suddenly drops for the largest box. In simulation 
set B (red), the quantity saturates for boxes with a radial extent 
that is twice as large as the azimuthal extent or larger. When 
keeping the radial extent constant (simulation set C, yellow line), 
the maximum saturates for the cubic box case. Hence, the only 
in azimuthal direction extended boxes do not lead to an artificial 
enhancement of the dust overdensities. The very high overdensity 
for run y-S seems to be a stochastic coincidence (compare top panel 
in Figure 0. The lower panel shows the dust overdensity against 
the correlation time. The measured points can be approximated 
with a power law (shown as a red solid line) with an exponent of 
0.38 ±0.05. The shaded region, with the dashed red lines as edges, 
gives the uncertainty of the fit. The one cross off the fit shows 
the results for run y-S. The blue square marks its position, if we 
neglect the two maxima shown in Figure It overlaps well with 
the fit region. 

weaker. Hence their radial drift velocity is almost not 
affected by the turbulence and they do not show strong 
concentrations. Small particles with low Stokes numbers 
are stronger coupled to the gas and, thus, also drift very 
slow. Particles with St ~ 1 are concentrated most by the 
zonal flow and, thus, have a stronger decreased radial 
velocity. Thus, the accumulation of dust particles is ex- 
pected to be strongest for particles with Stokes numbers 
around unity. For St = 0.01 particles, the drift velocity 
is strongly determined by the gas flow. This explains the 
strong deviation from the expected drift velocity. 

The upper right panel shows the total particles over- 
density normalized to the initial particle number den- 
sity. For run LspecMR (black diamonds), the smallest 
particles have higher concentrations than in the other 
simulations. This resulted from the choice of too few 
particles per grid cell. There only 100,000 particles per 
size bin were simulated. This results in overestimation, 
because the number density is normalized with the ini- 
tial number density no- For example, run LspecMRs (red 



squares) follows 2,000,000 particles per particle size bin. 
The highest concentrations were reached for particles of 
sizes St = 0.75 ... 5, as expected. However, the exact 
peak has a stochastic factor to it. Thus, the simulations 
peak at different particle sizes. The overdensities are 
more investigated in the lower row of panels. 

The surface number density of the particles is shown in 
the lower left panel. Here, the particles were integrated 
in the vertical direction. The trend is similar to the upper 
right panel. We read from this plot that particles with 
St = 0.1 are concentrated about ten times the initial con- 
centration. Together with the vertical overdensity due to 
sedimentation (lower right panel), a total overdensity of 
about 100 is created for St = 0.1 particles. 

The peaks in the vertical density structure of the par- 
ticles are shown in the lower right panel of Figure [121 
The Stokes number, St = Tff) defines the timescale after 
which the particles are settled down to the mid-plane. 
Particles with a high Stokes number are not fully settled 
down to the mid- plane, not even in the long-integration 
run MspecMRb. The resolution also limits this measure- 
ment for particles that are very close to the mid-plane. 
Smaller particles are not that strongly stratified. Thus, 
the vertical (Gaussian) structure is wider and shallower. 
This results in a lower value in this plot. The points for 
Stokes numbers 0.01-1 follow a power law with the index 
of 0.58 ± 0.03. The measured power law i ndex is slightly 
higher than the expected value of 0.5 ([Dubrulle et al.l 
119951 ). Most of the particles with St > 1 sediment very 
close to the mid-plane. This prohibits a further increase 
in the vertical density. A higher resolution and a mea- 
surement of the dust scale height is achieved in the next 
section. 

4.4.2. Dust Pressure Scale Height 

With a stratified particle distr ibution we can test th e 
vertical diffusion model (see, e.g. JCarballido et al.ll2006l ). 
The dust pressure scale height can be directly calcu- 
lated from the vertical positions of the particles of the 
same size. It is ap proximately proportion a l to S t~°' 5 
in agreement with Carballid o et al.l (|2006l l2011f ) and 
lYoudin fc Lithwickl ( 2007f ). The results are summarized 
in the upper panel of Figure HH Since the analyti- 
cal value was calculated with the a- value, the vertical 
Schmidt number 



Sc z = 



p, expected 



Hp, measured \Dt(oo) 



(16) 



can be calculated. We measured the vertical Schmidt 
number to have a very weak dependence on the particle 
size. In the lower panel of Figure [T3l we show that Sc z = 
3.4 -St - 11 . 

5. DISCUSSION AND CONCLUSIONS 

5.1. Zonal Flows and Axisymmetric Pressure Bumps 

Our simulations have dimensionless units. This allows 
us to interpret our results manyfold. We can pick the 
distance to the star in a certain range. In Section f275| we 
defined the global pressure gradient to be Av = 0.05c s . 
In the minimum mass solar nebula (MMSN) model, we 
can choose the distance t o the star to be between 0.35 
and 40 AU (Hayashi 1981). For this discussion, we pick 
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Figure 9. Top row shows the evolution of the azimuthal gas velocity and the dust density evolution of run L, respectively. The quantities 
are averaged in vertical and azimuthal direction and plotted in radial direction over time. The black dots in the upper left panel show the 
position of the highest dust density at each orbit. This shows that the overdensities of the dust often appear at places and times where the 
azimuthal gas velocity is high. This relation shows that the zonal flows accumulate dust and are a possible venue of planetesimal formation. 
The bottom left panel shows a scatter plot of the dust density p p against the azimuthal gas velocity, where both values are averaged in 
vertical and azimuthal direction, as in the upper panels. The bottom right panel shows a plot of the particle surface density p p (x,y,t) 
in relation to the azimuthal gas velocity, averaged in vertical direction, computed from a snapshot taken at 85T or b, the time when the 
maximum in the dust density occurs. Both plots show that it is more likely to find a high dust density at a location where the azimuthal 
gas velocity is high. 



r = 5 AU. In a thin disk model, we get a ratio for H/r — 
0.033(5 AU/AU) 1 / 4 -0.05; this defines us H = 0.25 AU. 
The isothermal sound speed is c s = HVt —66,000 cm s _1 . 
Thus, turbulent velocities (w rms ) are about 9000 cm s _1 
(—7000 cm s _1 for the high-resolution run LspecHR). 

Figure [14] shows the highest azimuthal velocity for all 
simulation sizes. We averaged over several maxima of 
u y (x, t) for every simulation to smooth over outliers. The 
zonal flows are super-Keplerian for all but runs XXL, x-S, 
and y-XL. In the largest box the flow only reaches slightly 
sub-Keplerian velocities. However, particles still get cap- 
tured in the resulting axisymmetric pressure bumps. The 
speeds measured in the largest simulation match those 



measured in lFlock et al.1 (|2011f ). 

We measured the radial size of the axisymmetric pres- 
sure bumps to be between 5 and 7H (see Figure [6]). At 
a distance of 5AU to the star, this size corresponds to 
— 1.25 ... 1.75 AU radial size for zonal flows, i.e., the dis- 
tance bet ween peaks o f (p) vz . This measurement agrees 
well with iSimon et al.l (|2012[ ) who measured the radial 
size of their zonal flows to be 6H. Further studies with 
varying box size in smaller steps could potentially narrow 
down the radial scale. 

We measured the lifetimes of the zonal flows up to 
50T or h. This agrees w ell with earlier stated lifetimes 
(jJohansen et al.ll2QQ9al : lUribe et al.ll2011f ). The strength 
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Figure 10. Radial drift velocity of the particles for two different 
simulations is shown in the upper panel. Particles in these simula- 
tions all have a Stokes number of St = Tff2 = 1. The orange lines 
are the exact measured radial velocities, averaged over all parti- 
cles. The black line represents the same value smoothed over the 
time of one local orbit. The blue dashed line shows the analytical 
resu lt in a stationary box for particles of St = 1 following Equation 
(|12|) . Particles in turbulent simulations generally drift slower than 
expected from the stationary solution, but the box size has little 
effect on the drift velocity. This is shown in the lower panel where 
the mean of the radial drift velocity is plotted against the box size. 
We omitted simulation set C for visibility. There is a minimum 
in drift speed for run L. However, this minimum is within the er- 
ror margins. The smallest errors are with run XXL; here the drift 
velocity drops by 28%. 

of the density bump reaches 15% and goes down to about 
10% in the largest simulation. The lower amplitude is 
consistent with the results from global simulations (pri- 
vate commun ication with Mario Floc k about the sim- 
ulations from iFlock et al.l I2011L 12012 ) who measured a 
density en hancement of slightly less than 10% . Some 
works (e.g. JUribe et al.ll2011c iSimon et al1l2012D measure 
stronger density enhancements. A possible explanation is 
that their a values are higher than in this work. Further 
studies on the dependence of volume average quantities 
to strength of zonal flows would be interesting. 

5.2. Dust in Zonal Flows 

Particles get trapped downstream of pressure bumps 
and build up over densities. To compare our dimension- 
less particle sizes with collision experiments and observa- 
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Figure 11. Distribution of the dust-to-gas ratio of the surface 
densities for runs S, M, L, and XL. For comparison a Poisson dis- 
tribution is shown with crosses. The initial distribution of the nu- 
merical simulation (dotted line) fits very well to the normal distri- 
bution. The average strength of the clustering for the dust surface 
density does not depend on the simulation box size. 

tions we have to assume a distance to the star and pick 
a solar system model. This will allow us to discuss our 
results in context to recent experiments. 

By choosing a model for the solar system, we can con- 
vert the dimensionless Stokes number St = nQ to a real 
particle size. The friction time n correlates to the parti- 
cle radius a with 



r (Ep) oy 

'f "^gas 



27rp. 



for Epstein drag and 



a = 



9r f (St) ^/i# 

4p,cr mo i 



(17) 



(18) 



for Stokes drag (see supplementary info for 
I Johansen et al.l l2007| ). Here S gas is the column 
density of the gas, p m is the density of solid material, 
fi = 3.9 x 10 _24 g is the mean molecular weight, 
and cr mo i = 2 x 10 _15 cm 2 i s the molecular cross 
section of molecular hydr ogen (jNakagawa et al.l 119861 : 
iChapman fc Cowlind[l97Ql ). 

The Epstein regime applies, if the particl e radius a 
does not exceed (9/4) (|Weidenschillinglll977af ) of the gas 
mean-free path 



A 



27T/J.H 
Pgas^mol ^gas^mol 



(19) 



The gas density and hence also the particle size for a 
given Stokes number St = Tffi depends very much on the 
used model. In Fig ure [T5l we overview four different mod- 
els. The MMSN (jWeidenschillinail977bl : [Havashil [l98lh 
was calculated from the mass of the existing planets, ne- 
glecting migration. Because this model allows no mass 
loss through accretion, often 3-MMSN is used to account 
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Figure 12. These panels show the behavior of particles with Stokes numbers from 0.01 to 100. The upper left panel shows the negative 
of the radial drift velocity and the relative difference between the measured and expected drift velocity. The dashed gray line shows the 
stationary solution for the radial drift, following Equation (1 1 2 f) . The highest drift velocities are obtained for particles with St = 1, but they 
are also slowed down strongest by the MRI-turbulence. The upper right panel shows the highest overdensity that occurred for the specific 
particle size during the entire simulation. The slopes for the different simulations match very well, apart from a jump around St = 1 (for 
LspecMRs and MspecMRb) and an offset for run LspecMR at small particle sizes. The former can be explained with the usage of a smaller 
simulation box for run MspecMRb (2.64 3 with weaker zonal flows) than in the other simulations (5.28 2 x 2.64 with stronger zonal flows). 
The offset showed that the number of particles per particle size was not sufficient in run LspecMR (10 5 particles in ~ 1.5 x 10 6 grid cells 
leads with five particles in one grid cell to a result of max (n p )/ (n p (t = 0)) = 75). The lower left panel shows the maximum of the column 
density for each particle size. It peaks at sizes of around St = 1. The lower right panel shows the maxima of the vertical distribution of 
particles. The cu rves (for St = 0.01 . . . 1) follow a power-law with the index of 0.58 zb 0.03. This is slightly steeper than the expected power 
law index of 0.5 (Dubrulle et al. 1995). In all four plots, the results of particles with St = 0.01 are to be interpreted with caution, because 
the simulations lacked sufficient amount of super-particles for these size bins. Further, large particles (St = 100) did not have enough time 
to sediment to the mid-plane. 



fo r some accretion. A low-density model was published 
bv lBrauer et al.l(|20Q8[ ). This model is adopted from mea- 
surements that indicate a s hallow surface densit y pro- 
file for protoplanetary disks (jAndrews et al. 2010). The 
high-density model was adopted from[Desch (2007), who 
introduced a "revised MMSN model" by using the start- 
ing positions in the Nice model of planetary dynamics 



([Tsiganis et al.l [20051 ). This model also takes planetary 
migration into account. The equations used to calculate 
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Figure 13. Dust scale height as a function of the Stokes number 
is shown in the upper panel. Different symbols depict the dif- 
fer entsjmulati^ dust scale height is calculated 
after Carballido et al. (2006) and compared with the fitted func- 
tion. This value is the vertical Schmidt number Sc z , shown in the 
lower panel. Its dependence on the particle size is very weak. 

the particle sizes in Figure [T5l are 
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(20) 



Throughout the discussion, we assume the MMSN model 
at 5 AU distance to the star for size reference for our test 
particles. This choice affects only the translation from 
the Stokes number St to a size, not the dynamics in our 
models. 

If the local dust density exceeds the Roche density, a 
clump is gravitationally bound against shea r. The Roche 
density can be approximated (|Kopallll989[ ) by 



PRoche(# = 5AU) = 



n 2 



4tt G(R 
100p(R = 5 AU) , 



5 AU) 



(21) 



for an MMSN. G is th e gravitational constant . 
The streaming instabilit y (jYoudin fc Johansenl 120071 : 
iJohansen fe Youdinl [2007) starts to act at dust-to-gas 
ratios of order unity. We started all our simulations 
with eo = p p (t = 0)/p = 0.01. Thus, a concentration 
of max(n p )/(no} = 100 corresponds to e s treaming = 1- 
The Roche density at 5 AU in an MMSN can be ex- 
pressed as eR OC he = PRoche/p ^ 100. We can see that 
objects of several decimeters up to some meters reach 
eRoche, while pebbles of some centimeters up to a decime- 
ter reach e str eaming from combining Figures [12] and [T5j 



0.10 

^ 0.08 
< 

+ 0.06 



x 



0.04 

0.02 
0.00 





I ] 


[ ; 


" ' ■ ] 


[ 


I ; 


I 




I " 




I 






I 


i ~ 









* 4r < -k ^Vv^VA 



Figure 14. Measured highest azimuthal velocity of all box size 
simulated. The gray dashed line shows the threshold to Keplerian 
velocity. Only runs XXL and y-XL do never get super-Keplerian. 
Simulation x-S does get super-Keplerian at some time, but not 
often enough to be significant. 



The concentration factors of run LspecHR in the upper 
right panel of Figure [T2l show us that with an initial dust- 
to-gas ratio of eo = 10 -2 particles of sizes St = 0.5 ... 10 
(St = 0.1 . . . 0.25) reach a dust-to-gas ratio of 100 (> 1). 
These sizes translate to 30 . . . 400 cm (6 ... 15 cm) in an 
MMSN at a 5 AU orbit using Figure [15l Considering 
back-reaction from the dust to the gas would allow the 
streaming instability to act. This will be subject of a 
future study. In our simulations, we see that the den- 
sity of 15 . . .600 cm sized icy boulders increases several 
thousand times over the equilibrium density, even with- 
out streaming instability and self-gravity of the particles. 
Sedimentation to the mid-plane leads to overdensities of 
around 40, while the contribution from the turbulence 
concentrates the boulders several hundred times. 

Since we do not study the influence of the back-reaction 
from particles to the gas, we were able to study several 
particle sizes in one simulation. That also means that 
the initial dust-to-gas ratio (eo) can be set arbitrary. We 
can interpret our results in the light of different metallic- 
ities. Particles with St > 0.5 will trigger the streaming 
instability even with eo = 10 -4 , while St = 0.1 particles 
need eo = 10 -2 . 

At the assumed distance in this discussion, the result- 
ing rings of trapped dust are not observable with cur- 
rent telescopes. If zonal flows form at larger distances 
to the star and dust rings form at an observable size, 
they could potentially be observable with ALMA. For an 
analysis one would have to adjust the parameter Av to 
account for the steeper pressure gradient. A preliminary 
study showed that particles of about 10 cm in size can 
get capture for a short amount of time at 100 AU dis- 
tance. However, this question goes beyond the scope of 
this paper and should be addressed in a future study. 

6. SUMMARY AND OUTLOOK 

We performed numerical simulations of MRI-driven 
turbulence in shearing boxes, covering the parameter 
space for radial and azimuthal box sizes up to 21.12i^. 
Further, we followed the reaction of the dust particle den- 
sity to the turbulence. Our major findings are as follows. 

1 . Turbulent energy and stresses double when increas- 
ing the azimuthal size of the simulation from 1.32 
to 2.64 pressure scale heights. Turbulence pa- 
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Figure 15. Particle sizes as a function of the dimensionless Stokes 
number for the four discussed models at 5 AU. The black squares 
show the used Stokes numbers and their corresponding size in the 
case of the minimum mass solar nebula model. 



rameters in radially small box sizes stay approx- 
imately constant. This confirms the results in 
Fr omang fc Stond ([20091 ). In larger boxes, turbu- 
lent fluctuations and stresses are observed to re- 
main ^^onstajr^a^amst^l^ges in the box size (see 
also Uohansen et al.l 12009a). This rap i d con ver- 
gence was also observed in iSimon et al.l (|2012[ ) . 



2. Surface density fluctuations grow to large scales 
in the box and have lifetimes of up to 50 or- 
bits. The scales of these pressure bumps increase 
with increasing radial box size, until it saturates 
at approximately 5-7 pressure scale heights. The 
scales are decreased when the azimuthal box size 
is much more increased than the radial box size. 
The radial scales of the pressure bumps are con- 
sistent with the len gth sc ales measured in local 
(e.g.. Uohansen et al.ll 2009a: Simon et al. 2012) and 
global (e.g., iLvra et al.l 120081 : lUribe et al.l 120111 1 
simulations. This might be the natural size of 
these over densities. The pressure bumps are in 
geostrophic balance with sub- and super-Keplerian 
zonal flows. At 5 AU distance to the star 6H corre- 
spond to ~ 1.5 AU. The amplitude of the density 
bump reaches 15% and goes down to about 10% in 
the largest simulation. 

3. Particles with St = Tffl = 1 are getting trapped 
efficiently by the axisymmetric pressure bumps. 
They accumulate in regions of minima in the sec- 
ond derivative of the gas density as predicted ana- 
lytically (e.g. JKlahr fcT53l2QQlh . The concentra- 
tion factor correlates with the correlation time of 
the zonal flows. Hence, the first two steps of plan- 
etesimal formation^ in protoplanetary disk with an 
acting MRI are: vertical settling via sedimentation 

7 After coagulation from /im-sized particles to St = 0.1, 1. 



and radial concentration by trapping of dust in ax- 
isymmetric pressure bumps. Further concentration 
comes likely from stochastic processes. Clustering 
properties do not depend strongly on strength or 
lifetime of the zonal flows. 

4. We reach dust-to-gas ratios of 50-100. These den- 
sities are of the order of the Roche density at 5 AU 
in an MMSN. The dust overdensities scale with the 
lifetime of the zonal flow structures by a power law 
with an exponent of 0.38 ± 0.05 (see Figure [8]). To 
what degree these high dust-to-gas ratios disturb 
the axisymmetric pressure bumps that developed 
in the zonal flows has to be investigated in further 
studies with back-reaction to the gas. 

5. Particles of only a few centimeters in size (at 5 AU 
in an MMSN, St = 0.1) accumulate in overdensi- 
ties that are increased by a factor of ~ 100, lead- 
ing to a dust-to-gas ratio of 1 in the mid-plane, 
thus triggering the streaming instability. With- 
out MRI and zonal flows St = 0.1 particles do 
not clump strongly and cannot trigger the stream- 
ing instability for solar metallicity Z = eo = 0.01 
(jJohansen et al.ll2009b[ ). 

This is the first work on the effect from large-scale 
zonal flows on dust particles in an MHD simulation. Dust 
gets trapped downstream of long-lived high-pressure re- 
gions and achieves overdensities that have the potential 
to generate streaming instability and to become gravita- 
tionally unstable. Planetesimal formation in large boxes 
will be further investigated in simulations with particle 
feedback on the gas and self-gravitating particles in a 
future study. 

In the future, we will focus on one model and study 
various initial dust-to-gas ratios and particle size distri- 
butions. We will probably use the already converged 
run L (5.28H x 5.2SH x 2.64 Jf). This choice is also a 
trade-off between simulation box size and computational 
expense. 

K.D. and H.K. have been supported by the Deutsche 
Forschungsgemeinschaft Schwerpunktprogramm (DFG 
SPP) 1385 "The first ten million years of the solar sys- 
tem"; K.D. received further support by the IMPRS for 
Astronomy & Cosmic Physics at the University of Hei- 
delberg. A.J. was partially funded by the European 
Research Council under ERC Starting Grant agreement 
278675PEBBLE2PLANET and by the Swedish Research 
Council (grant 20103710). Computer simulations have 
been performed at the THEO cluster at Rechenzentrum 
Garching and at the JUGENE cluster in Forschungszen- 
trum Julich (project number HHD19). We thank an 
anonymous referee for a very thorough referee report that 
improved the quality of the paper greatly. 

REFERENCES 



Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C, & 

Dullemond, C. P. 2010, ApJ, 723, 1241 
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 
Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1 
Beitz, E., Guttler, C, Blum, J., et al. 2011, ApJ, 736, 34 
Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 
Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 



18 



Dittrich, Klahr, & Johansen 



Brandenburg, A., Nordhmd, A., Stein, R. F., & Torkelsson, U. 

1995, ApJ, 446, 741 
Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 
Carballido, A., Bai, X.-N., & Cuzzi, J. N. 2011, MNRAS, 415, 93 
Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 

1633 

Cassen, P., & Moosman, A. 1981, Icar, 48, 353 

Chapman, S., & Cowling, T. G. 1970, The Mathematical Theory 
of Non-uniform Gases. An Account of the Kinetic Theory of 
Viscosity, Thermal Conduction and Diffusion in gases (3rd ed.; 
Cambridge: Univ. Press) 

Cuzzi, J. N., Hogan, R. C, & Shariff, K. 2008, ApJ, 687, 1432 

Desch, S. J. 2007, ApJ, 671, 878 

Dominik, C, Blum, J., Cuzzi, J. N., & Wurm, G. 2007 in 

Protostars and Planets V ed. B. Reipurth, D. Jewitt, &; K. Keil 

(Tucson, AZ: Univ. Arizona Press), 783 
Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icar, 114, 237 
Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, 

T. 2010, A&A, 515, A70 
Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, 

R., & Oliveira, J. M. 2010, A&A, 510, A72 
Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., &; Henning, 

T. 2011, ApJ, 735, 122 
Flock, M., Dzyurkevich, N., Klahr, H., Turner, N., & Henning, T. 

2012, ApJ, 744, 144 
Fromang, S., & Stone, J. M. 2009, A&A, 507, 19 
Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, 

ApJ, 694, 1010 

Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJL, 553, 
153 

Hayashi, C. 1981, PThPS, 70, 35 

Heinemann, T., & Papaloizou, J. C. B. 2009, MNRAS, 397, 64 
Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292 
Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62 
Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Natur, 
448, 1022 

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 
Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269 
Johansen, A., Youdin, A., & Mac Low, M.-M. 2009b, ApJL, 704, 
75 

Klahr, H. H., & Lin, D. N. C. 2001, ApJ, 554, 1095 



Kopal, Z. 1989, The Roche Problem and its Significance for 
Double-star Astronomy (Astrophysics and Space Science 
Library, Vol. 152; Dordrecht: Kluwer) 

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71 

Lodders, K. 2003, ApJ, 591, 1220 

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 
479, 883 

Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icar, 67, 375 
Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G., & Norman, M. L. 

2011, ApJ, 740, 6 

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 
Safronov, V. S. (ed.) 1969, Evoliutsiia doplanetnogo oblaka 

(English transl.: Evolution of the Protoplanetary Cloud and 

Formation of Earth and the Planets, NASA Tech. Transl. 

F-677, 1972,; Jerusalem: Israel Sci. Transl.) 
Shakura, N. L, & Sunyaev, R. A. 1973, A&A, 24, 337 
Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 

422, 2685 

Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 

2012, ApJ, 749, 189 

Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142 
Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, 
Natur, 435, 459 

Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 
85 

von Neumann, J., & Richtmyer, R. D. 1950, JAP, 21, 232 
Weidenschilling, S. J. 1977a, MNRAS, 180, 57 
Weidenschilling, S. J. 1977b, Ap&SS, 51, 153 
Weidenschilling, S. J. 1997, Icar, 127, 290 

Whipple, F. L. 1972, in Proc. Twenty-First Nobel Symp. on From 
Plasma to Planet, ed. A. Evlius (New York: Wiley), 211 

Windmark, F., Birnstiel, T., Guttler, C, et al. 2012a, A&A, 540, 
A73 

Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 

2012b, A&A, 544, L16 
Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 
Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 
Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588 
Zsom, A., Ormel, C. W., Guttler, C, Blum, J., & Dullemond, 

C. P. 2010, A&A, 513, A57 



